home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_02 / qsieve.c < prev    next >
Text File  |  1989-04-17  |  13KB  |  465 lines
  1. /*
  2.  *   Program to factor big numbers using Pomerance-Silverman-Montgomery  
  3.  *   multiple polynomial quadratic sieve.
  4.  *   See "The Multiple Polynomial Quadratic Sieve", R.D. Silverman,
  5.  *   Math. Comp. Vol. 48, 177, Jan. 1987, pp329-339
  6.  *
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include <math.h> 
  11. #include "miracl.h"
  12.  
  13. #define MEM 240        /* Maximum size of factor base   */
  14. #define MLF 100        /* Max. number of large factors  */
  15. #define PAK 1+MEM/(2*BTS)
  16. #define SSIZE 10000    /* Maximum sieve size            */
  17. #define LOG log        /* if log not available use any other log */
  18.  
  19. big NN,TT,DD,RR,VV,PP,XX,YY,DG,IG,AA,BB;
  20. big x[MEM+1],y[MEM+1],z[MLF+1],w[MLF+1];
  21. int ep[MEM+1],r1[MEM+1],r2[MEM+1],rp[MEM+1],pr[MLF+1];
  22. int b[MEM+1],e[MEM+1],hash[2*MLF+1];
  23. unsigned int EE[MEM+1][PAK],G[MLF+1][PAK];
  24. unsigned char logp[MEM+1],sieve[SSIZE+1];
  25. int mm,jj,nbts,nlp,lp,hmod,hmod2;
  26. bool partial;
  27.  
  28. extern char *malloc();
  29.  
  30. int knuth(ep,N,R,D,T)
  31. int ep[];
  32. big N,R,D,T;
  33. { /* initialization */
  34.     double fks,dp,top;
  35.     bool found;
  36.     int *primes;
  37.     int i,j,mm,bk,nk,m,d,kk,r,p;
  38.     static int K[]={0,1,2,3,5,6,7,10,11,13,14,15,17,0};
  39.     printf("input number to be factored N= \n");
  40.     d=cinnum(N,stdin);
  41.     if (prime(N))
  42.     {
  43.         printf("this number is prime!\n");
  44.         exit(0);
  45.     }
  46.     if (d<10) m=d;
  47.     else m=25;
  48.     if (d>20) m=30*(d-20);
  49.     if (m>MEM) m=MEM;
  50.     primes=gprime((-2)*(m+5));
  51.     top=(-10.0e0);
  52.     found=FALSE;
  53.     nk=0;
  54.     bk=0;
  55.     ep[0]=1;
  56.     ep[1]=2;
  57.     do
  58.     { /* search for best Knuth-Schroepel multiplier */
  59.         mm=m;
  60.         kk=K[++nk];
  61.         if (kk==0)
  62.         { /* finished */
  63.             kk=K[bk];
  64.             found=TRUE;
  65.         }
  66.         premult(N,kk,D);
  67.         if (root(D,2,R))
  68.         {
  69.             printf("%dN is a perfect square!\n",kk);
  70.             continue;
  71.         }
  72.         fks=LOG(2.0e0)/(2.0e0);
  73.         r=subdiv(D,8,T);
  74.         if (r==1) fks*=(4.0e0);
  75.         if (r==5) fks*=(2.0e0);
  76.         fks-=LOG((double)kk)/(2.0e0);
  77.         i=0;
  78.         j=1;
  79.         while (j<m)
  80.         { /* select small primes */
  81.             p=primes[++i];
  82.             if (p==0)
  83.             {
  84.                 mm=j;
  85.                 break;
  86.             }
  87.             r=subdiv(D,p,T);
  88.             if (spmd(r,(p-1)/2,p)<=1)
  89.             { /* use only if Jacobi symbol = 0 or 1 */
  90.                 ep[++j]=p;
  91.                 dp=(double)p;
  92.                 if (kk%p==0) fks+=LOG(dp)/dp;
  93.                 else         fks+=2*LOG(dp)/(dp-1.0e0);
  94.             }
  95.         }
  96.         if (fks>top)
  97.         { /* find biggest fks */
  98.             top=fks;
  99.             bk=nk;
  100.         }
  101.     } while (!found);
  102.     printf("using multiplier k= %d\n",kk);
  103.     printf("and %d small primes as factor base\n",mm);
  104.     free(primes);
  105.     return mm;
  106. }
  107.  
  108. bool factored(lptr,T)
  109. big T;
  110. long lptr;
  111. { /* factor quadratic residue */
  112.     bool facted;
  113.     int i,j,r,st;  
  114.     partial=FALSE;
  115.     facted=FALSE;
  116.     for (j=1;j<=mm;j++)
  117.     { /* now attempt complete factorisation of T */
  118.         r=lptr%ep[j];
  119.         if (r<0) r+=ep[j];
  120.         if (r!=r1[j] && r!=r2[j]) continue;
  121.         while (subdiv(T,ep[j],XX)==0)
  122.         { /* cast out ep[j] */
  123.             e[j]++;
  124.             copy(XX,T);
  125.         }
  126.         st=size(T);
  127.         if (st==1)
  128.         {
  129.            facted=TRUE;
  130.            break;
  131.         }
  132.         if (size(XX)<=ep[j])
  133.         { /* st is prime < ep[mm]^2 */
  134.             if (st>=TOOBIG || (st/ep[mm])>=(1+MLF/50)) break;
  135.             if (st<=ep[mm])
  136.                 for (i=j;i<=mm;i++)
  137.                 if (st==ep[i])
  138.                 {
  139.                     e[i]++;
  140.                     facted=TRUE;
  141.                     break;
  142.                 }
  143.             if (facted) break;
  144.             lp=st;  /* factored with large prime */
  145.             partial=TRUE;
  146.             facted=TRUE;
  147.             break;
  148.         }
  149.     }
  150.     return facted;
  151. }
  152.  
  153. bool gotcha()
  154. { /* use new factorisation */
  155.     int r,j,i,k,n,rb,had,hp;
  156.     bool found;
  157.     found=TRUE;
  158.     if (partial)
  159.     { /* check partial factorisation for usefulness */
  160.         had=lp%hmod;
  161.         forever
  162.         { /* hash search for matching large prime */
  163.             hp=hash[had];
  164.             if (hp<0)
  165.             { /* failed to find match */
  166.                 found=FALSE;
  167.                 break;
  168.             }
  169.             if (pr[hp]==lp) break; /* hash hit! */
  170.             had=(had+(hmod2-lp%hmod2))%hmod;
  171.         }
  172.         if (!found && nlp>=MLF) return FALSE;
  173.     }
  174.     copy(PP,XX);
  175.     convert(1,YY);
  176.     for (k=1;k<=mm;k++)
  177.     { /* build up square part in YY  *
  178.        * reducing e[k] to 0s and 1s */
  179.         if (e[k]<2) continue;
  180.         r=e[k]/2;
  181.         e[k]%=2;
  182.         convert(ep[k],TT);
  183.         power(TT,r,TT,TT);
  184.         multiply(TT,YY,YY);
  185.     }
  186. /* debug only - print out details of each factorization 
  187.     cotnum(XX,stdout);
  188.     cotnum(YY,stdout);
  189.     if (e[0]==1) printf("-1");
  190.     else printf("1");
  191.     for (k=1;k<=mm;k++)
  192.     {
  193.         if (e[k]==0) continue;
  194.         printf(".%d",ep[k]);
  195.     }
  196.     if (partial) printf(".%d\n",lp);
  197.     else printf("\n");
  198. */
  199.     if (partial)
  200.     { /* factored with large prime */
  201.         if (!found)
  202.         { /* store new partial factorization */
  203.             hash[had]=nlp;
  204.             pr[nlp]=lp;
  205.             copy(XX,z[nlp]);
  206.             copy(YY,w[nlp]);
  207.             for (n=0,rb=0,j=0;j<=mm;j++)
  208.             {
  209.                 G[nlp][n]|=((e[j]&1)<<rb);
  210.                 if (++rb==nbts) n++,rb=0;
  211.             }
  212.             nlp++;
  213.         }
  214.         if (found)
  215.         { /* match found so use as factorization */
  216.             printf("\b\b\b\b\b*");
  217.             mad(XX,z[hp],XX,NN,NN,XX);
  218.             mad(YY,w[hp],YY,NN,NN,YY);
  219.             for (n=0,rb=0,j=0;j<=mm;j++)
  220.             {
  221.                 e[j]+=((G[hp][n]>>rb)&1);
  222.                 if (e[j]==2)
  223.                 {
  224.                     premult(YY,ep[j],YY);
  225.                     divide(YY,NN,NN);
  226.                     e[j]=0;
  227.                 }
  228.                 if (++rb==nbts) n++,rb=0;
  229.             }
  230.             premult(YY,lp,YY);
  231.             divide(YY,NN,NN);
  232.         }
  233.     }
  234.     else printf("\b\b\b\b\b ");
  235.     if (found)
  236.     {
  237.         for (k=mm;k>=0;k--)
  238.         { /* use new factorization in search for solution */
  239.             if (e[k]%2==0) continue;
  240.             if (b[k]<0)
  241.             { /* no solution this time */
  242.                 found=FALSE;
  243.                 break;
  244.             }
  245.             i=b[k];
  246.             mad(XX,x[i],XX,NN,NN,XX);    /* This is very inefficient -  */
  247.             mad(YY,y[i],YY,NN,NN,YY);    /* There must be a better way! */
  248.             for (n=0,rb=0,j=0;j<=mm;j++)
  249.             { /* Gaussian elimination */
  250.                 e[j]+=((EE[i][n]>>rb)&1);
  251.                 if (++rb==nbts) n++,rb=0;
  252.             }
  253.         }
  254.         for (j=0;j<=mm;j++)
  255.         { /* update YY */
  256.             if (e[j]<2) continue;
  257.             convert(ep[j],TT);
  258.             power(TT,e[j]/2,NN,TT);
  259.             mad(YY,TT,YY,NN,NN,YY);
  260.         }
  261.         if (!found)
  262.         { /* store details in EE, x and y for later */
  263.             b[k]=jj;
  264.             copy(XX,x[jj]);
  265.             copy(YY,y[jj]);
  266.             for (n=0,rb=0,j=0;j<=mm;j++)
  267.             {
  268.                 EE[jj][n]|=((e[j]&1)<<rb);
  269.                 if (++rb==nbts) n++,rb=0;
  270.             }
  271.             jj++;
  272.             printf("%4d",jj);
  273.         }
  274.     }
  275.     if (found)
  276.     { /* check for false alarm */
  277.         printf("\ntrying...\n");
  278.         add(XX,YY,TT);
  279.         if (compare(XX,YY)==0 || compare(TT,NN)==0) found=FALSE;
  280.         if (!found) printf("working... %4d",jj);
  281.     }
  282.     return found;
  283. }
  284.  
  285. void initv()
  286. { /* initialize big numbers and arrays */
  287.     int i,j;
  288.     NN=mirvar(0); 
  289.     TT=mirvar(0);
  290.     DD=mirvar(0);
  291.     RR=mirvar(0);
  292.     VV=mirvar(0);
  293.     PP=mirvar(0);
  294.     XX=mirvar(0);
  295.     YY=mirvar(0);
  296.     DG=mirvar(0);
  297.     IG=mirvar(0);
  298.     AA=mirvar(0);
  299.     BB=mirvar(0);
  300.     for (i=0;i<=MEM;i++)
  301.     {
  302.         x[i]=mirvar(0);
  303.         y[i]=mirvar(0);
  304.     }
  305.     for (i=0